---
title: "Lab 11 Code"
author: "Fabio Setti"
date: "Spring-2025"
output: html_document
---


## Today's Packages and Data 🤗

```{r}
# install.packages("tidyverse")
install.packages("emmeans")
```

```{r}
library(tidyverse)
theme_set(theme_classic(base_size = 14, 
                        base_family = 'serif'))
library(emmeans)
```

```{r}
hsb2 <- rio::import("https://github.com/quinix45/PSYC-7804-Regression-Lab-Slides/raw/refs/heads/main/Slides%20Files/Data/hsb2.csv")
```


## What are we looking at Today?


```{r}
hsb2$socst_cent <- scale(hsb2$socst, scale = FALSE)[,1]

hsb2$gender <- factor(hsb2$gender)

contrasts(hsb2$gender)
```


## Starting from The Top


```{r}
reg_gender <- lm(write ~ gender, data = hsb2)

summary(reg_gender)
```



## Visualizing `write ~ gender`


```{r}
mean_female <- mean(hsb2$write[hsb2$gender == "female"])
mean_male <- mean(hsb2$write[hsb2$gender == "male"])

ggplot(hsb2, aes(x = gender, y = write)) +
  geom_point() +
  geom_hline(aes(yintercept = mean(mean_female), color = "Female"), 
             linetype = "dashed") +
  geom_hline(aes(yintercept = mean(mean_male), color = "Male"), 
             linetype = "dashed") +
  scale_color_manual(values = c("Female" = "blue", "Male" = "red")) +
  labs(color = "Means")
```



## Adding a continous Predictor


```{r}
reg_both <- lm(write ~ gender + socst_cent, data = hsb2)

summary(reg_both)
```



## Visualizing `write ~ gender + socst`


```{r}
ggplot(hsb2, aes(x = socst_cent, y = write, colour = gender)) +
  geom_point(aes(shape = gender), alpha= .6) +
  geom_abline(intercept = coef(reg_both)[1], slope = coef(reg_both)[3], col = "purple") + 
  geom_abline(intercept = coef(reg_both)[1] + coef(reg_both)[2], slope = coef(reg_both)[3], col = "red") + 
  scale_color_manual(values=c( "purple", "red")) 
```


## Adding the Interaction

```{r}
reg_int <- lm(write ~ gender*socst_cent, data = hsb2)

summary(reg_int)
```


## Visualizing `write ~ gender * socst`


```{r}
ggplot(hsb2, aes(x = socst_cent, y = write, colour = gender)) +
  geom_point(aes(shape = gender), alpha= .6) +
  geom_abline(intercept = coef(reg_int)[1], slope = coef(reg_int)[3], col = "purple") + 
  geom_abline(intercept = coef(reg_int)[1] + coef(reg_int)[2], slope = coef(reg_int)[3] + coef(reg_int)[4], col = "red") + 
  scale_color_manual(values=c( "purple", "red")) 
```


##  Simple Slopes With `emmeans`

```{r}
emtrends(reg_int, ~gender, var = "socst_cent")
```

## Predicted Means 

```{r}
sd_below <- -sd(hsb2$socst_cent)
sd_above <- sd(hsb2$socst_cent)

# the mean is 0 because the variable is centered
mylist <- list(socst_cent=c(sd_below,
                            0,
                            sd_above))
```


```{r}
means <- emmeans(reg_int, ~socst_cent*gender, 
                at = mylist)

means
```

## Predicted Means Graphically

```{r}
ggplot(hsb2, aes(x = socst_cent, y = write, colour = gender)) +
  geom_point(aes(shape = gender), alpha= .6) +
  geom_abline(intercept = coef(reg_int)[1], slope = coef(reg_int)[3], col = "purple") + 
  geom_abline(intercept = coef(reg_int)[1] + coef(reg_int)[2], slope = coef(reg_int)[3] + coef(reg_int)[4], col = "red") + 
  scale_color_manual(values=c( "purple", "red")) +
  geom_point(aes(x = -10.7, y = 50.3), color = "brown", size = 4) +
  geom_point(aes(x = -10.7, y =  43.8), color = "brown", size = 4) +
  geom_point(aes(x = 0, y = 54.8), color = "brown", size = 4) +
  geom_point(aes(x = 0, y =  50.5), color = "brown", size = 4) +
  geom_point(aes(x = 10.7, y = 59.3), color = "brown", size = 4) +
  geom_point(aes(x = 10.7, y = 57.2), color = "brown", size = 4) 
```


##  Testing Mean Differences

```{r}
contrast(means, "pairwise", by= "socst_cent")
```


## Mean differences Graphically


```{r}
sum <- summary(means)

ggplot(hsb2, aes(x = socst_cent, y = write, colour = gender)) +
  geom_point(aes(shape = gender), alpha= .2) +
  geom_abline(intercept = coef(reg_int)[1], slope = coef(reg_int)[3], col = "purple") + 
  geom_abline(intercept = coef(reg_int)[1] + coef(reg_int)[2], slope = coef(reg_int)[3] + coef(reg_int)[4], col = "red") + 
  scale_color_manual(values=c( "purple", "red")) +
  geom_segment(y = sum[1,3], yend = sum[4,3], x = sum[4,1], xend = sum[4,1], col = "black", lty = 2) +
  geom_segment(y = sum[2,3], yend = sum[5,3], x = sum[5,1], xend = sum[5,1], col = "black", lty = 2) +
  geom_segment(y = sum[6,3], yend = sum[3,3], x = sum[3,1], xend = sum[3,1], col = "black", lty = 2) 
  
```


## johnson-neyman Plot With Categorical Moderator?


```{r}
# the predictor must be a vector of 0s and 1s
hsb2$gender_bin <- ifelse( hsb2$gender == "female", 0, 1)

# rerun the regression with binary gender
reg_int_2 <- lm(write ~ gender_bin * socst_cent, data = hsb2)

# save summary for mean differences
interactions::johnson_neyman(reg_int_2, 
                             pred = "gender_bin",
                             modx = "socst_cent")
```



## Plotting with `emmip()`

```{r}
min <- min(hsb2$socst_cent)
max <- max(hsb2$socst_cent)

mylist <- list(socst_cent = c(min, 
                              max))


emmip(reg_int, gender~socst_cent, at = mylist)
```


## Plotting with `interact_plot()`

```{r}
interactions::interact_plot(reg_int_2, 
              modx = "gender_bin",
              pred = "socst_cent",
              modx.values = c(0, 1))
```



